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Abstract 


Molecular  dynamics  simulations  at  constant  temperature  have  been  carried 
out  for  the  primitive  model  of  1-3  electrolyte  solutions.  Thermodynamics,  pair 
distribution  functions,  and  self-diffusion  coefficients  were  computed  to  examine 
the  electrostatic  effects  on  the  structural  and  dynamical  properties.  The  simula¬ 
tion  results  were  used  to  evaluate  various  theoretical  equations,  namely,  the 
exponential  form  of  Debye-Huckel  theory,  the  mean  spherical  approximation,  and 
the  hypernetted  chain  approximation.  As  has  been  observed  for  symmetrical 
electrolytes,  the  latter  turns  out  to  be  the  best  approximation.  For  asymmetri¬ 
cally  charged  1-3  electrolytes,  it  was  found  that  ionic  aggregation  significantly 
influenced  the  dynamical  properties  of  electrolytes.  Coherent  motion  between 
highly  charged  negative  ions  and  positive  ions  surrounding  them  was  deduced 
from  the  time  dependence  of  the  velocity  autocorrelation  functions,  particularly 
at  concentrations  between  0.4  and  4  total  molar. 
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1.  Introduction 

During  the  last  couple  of  decades  much  progress  has  been  made  in  our  under¬ 
standing  and  interpretation  of  the  physical  properties  of  electrolyte  solutions.  A 
number  of  different  theoretical  approximations  [1-7]  have  been  used  to  evaluate 
the  thermodynamic  and  structural  properties  of  bulk  electrolytes.  One  of  the 
simplest  but  most  commonly  used  models  in  these  theoretical  approaches  is  the 
so-called  "primitive  model  electrolyte  solution"  [1.2],  in  which  the  charged  hard- 
sphere  ions  are  immersed  in  a  continuum  solvent  represented  only  as  a  uniform 
dielectric  constant  of  the  medium. 

Reliable  and  unambiguous  results,  in  turn,  have  become  increasingly  neces¬ 
sary  to  eliminate  any  underlying  uncertainties  involved  in  theoretical  approxima¬ 
tions.  However,  the  present  level  of  modeling,  particularly  the  use  of  a  solvent 
continuum  assumption,  is  far  too  crude  to  allow  the  direct  comparison  with  real 
laboratory  experiment.  Consequently,  machine  experiments  (computer  simula¬ 
tions),  which  provide  essentially  "exact"  experimental  data  for  precisely  defined 
model  systems,  have  proven  to  be  an  extremely  useful  diagnostic  tool  to  investi¬ 
gate  such  systems.  The  great  advantage  of  computer  simulations  over  experiment 
lies  in  the  possibility  of  obtaining  detailed  physical  information,  which  may  be 
very  difficult  or  impossible  to  realize  in  the  laboratory.  The  best  known  example 
is  the  measurement  of  pair  distribution  functions  as  a  function  of  distance  at  the 
molecular  level. 

There  are  in  general  two  classes  of  computer  experiment:  stochastic  Monte 
Carlo  (MC)  and  deterministic  molecular  dynamics  (MD)  simulations.  In  the  MD 
calculations,  the  actual  trajectories  of  molecules  are  evaluated  by  the  numerical 
integration  of  Newton’s  equations  of  motion.  In  addition  to  static  equilibrium 
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properties,  time-dependent  transport  properties  can  be  also  obtained  by  the  NED 
method.  For  the  system  of  primitive  model  electrolytes,  however,  the  discontinu¬ 
ous  nature  of  hard-sphere  repulsions  combined  with  continuous  electrostatic 
interactions  introduces  some  technical  difficulties  into  the  traditional  \1D  algo¬ 
rithm.  Recently,  implementations  have  been  made  to  extend  the  MD  method  to 
the  system  of  hard-cores  with  soft  attractive  potentials.  In  the  work  of  Heyes  [8] 
for  the  restricted  primitive  model  of  1:1  electrolytes,  two  different  forms  of  MD 
were  combined.  The  hard-sphere  collisions  were  allowed  to  take  place  while  the 
forces  and  velocity  changes  due  to  continuous  electrostatic  interactions  were 
evaluated.  The  agreement  with  the  previous  MC  calculations  was  excellent  even 
at  higher  concentrations. 

Almost  all  simulations  for  the  primitive  electrolytes  have  been  carried  out 
using  the  MC  simulation  method  [9-12].  Much  of  the  earlier  work  on  the  MC 
calculations  in  a  variety  of  different  ensembles  has  been  summarized  by  Levesque 
et  al.  [13].  Valleau  and  his  collaborators  have  extensively  reported  the  canonical 
and  the  grand  canonical  MC  results  [9]  of  such  systems.  In  their  studies,  various 
theoretical  approaches  were  also  discussed  and  compared  with  the  results 
obtained  from  their  MC  simulations.  They  adopted  the  minimum  image  conven¬ 
tion  to  evaluate  the  Coulombic  part  of  potentials  for  each  configuration,  and  the 
resulting  energy  was  linearly  extrapolated  as  a  function  of  l/N  to  estimate  the 
values  of  an  infinite  system.  These  data  should  be  accepted  with  some  care 
because  the  method  with  as  few  as  200  ions  the  method  could  yield  results 
depending  on  system  size.  This  is  particularly  true  for  higher  concentrations  and 
for  higher  charged  systems.  Explicit  results  for  osmotic  pressure  were  not  pub¬ 
lished  in  the  previous  works. 

In  the  MD  simulations,  better  statistics  can  be  achieved.  For  instance,  in 
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order  to  calculate  the  virial  contribution  to  the  equation  of  state  for  the  hard¬ 
core  system,  the  MC  computations  require  an  accurate  estimation  of  pair  distri¬ 
bution  functions' at  the  contact  point  [14].  For  the  system  of  ionic  solutions,  pair 
functions  may  change  rapidly  near  contact  distances  due  to  the  formation  of 
ionic  aggregates  (see  Fig.  4a).  The  extrapolation  of  pair  functions  to  the  contact 
value  lead  to  large  uncertainties  in  this  case  of  ionic  solutions.  For  this  reason  the 
MC  results  for  osmotic  pressure  coefficients  are  known  to  be  less  certain  than 
those  for  energy  calculations. 

In  the  present  paper,  we  report  computer  simulation  results  for  the  asym¬ 
metrically  charged  system  of  1-3  electrolyte  solutions  via  molecular  dynamics 
simulation  techniques.  This  asymmetric  system  is  of  special  interest  because  it 
provides  a  strong  test  of  the  applicability  of  approximate  theories  available  in  the 
literature.  Computational  details  of  the  method  employed  here  are  discussed  in 
the  next  section.  In  section  3,  the  thermodynamic  and  structural  properties 
obtained  from  our  MD  simulations  are  compared  with  various  theoretical  predic¬ 
tions,  namely,  the  exponential  form  of  Debye-Huckel  theory  [9,15],  the  mean 
spherical  approximation  [3],  and  the  hypernetted  chain  approximation  [4],  We 
also  report  in  this  section  the  dynamical  properties  including  the  velocity  auto¬ 
correlation  functions  and  the  self- diffusion  coefficients.  The  dynamical  properties 
of  the  primitive  model  electrolytes  are  of  course  not  comparable  to  those  of  elec¬ 
trolyte  solutions  since  solvent  dynamics  are  totally  absent  in  this  model.  How¬ 
ever,  the  results  are  an  inexpensive  by-product  of  the  simulation  and  should  be 
useful  for  testing  the  approximate  theories  of  transport  in  charged  systems  (e.g., 
plasmas  and  molten  salts). 


2.  Model  and  Computations 


The  MD  calculations  were  carried  out  with  a  system  containing  216  charged 
hard-spheres  with  equal  diameter,  d.  of  0.425  nm  and  mass,  m,  of  100  a.m.u. 
(There  are  162  cations  with  charge  +1  and  54  anions  with  charge  -3.)  Usual 
periodic  boundary  conditions  were  applied  in  a  cubic  fundamental  cell  to  approx¬ 
imate  an  infinite  system. 


The  pair  potential  energy  between  ions  i  and  j  is  (in  S.I.  units) 


tT„(r) 


DC 


if  x<d 


ZjZi<T 

47re0err 


if  x>d, 


(1) 


where  z,  and  z]  are  the  valence  of  ion  i  and  j,  and  q  is  the  unit  of  electronic 
charge.  eT  is  the  uniform  dielectric  constant  of  the  medium  relative  to  the  per¬ 
mittivity  of  free  space,  e0.  As  used  in  most  other  studies,  the  relative  dielectric 
constant,  er,  was  chosen  to  be  78.356  corresponding  to  those  for  water  at  room 
temperature,  298.16  K. 


The  long-ranged  interactions  in  ionic  systems  give  a  configurational  energy 
that  converges  slowly  with  the  increasing  system  size.  Use  of  a  spherical  cut-off 
or  a  minimum  image  convention  has  shown  to  be  inappropriate,  particularly  for 
highly  charged  dense  systems  [16].  The  resulting  configurational  energy  for  such 
systems  must  take  account  for  the  ion  pairs  not  only  with  the  nearest  images  in  a 
fundamental  cell  but  also  with  all  images  in  other  periodic  cells.  The  Coulombic 
potentials  or  forces,  in  this  study,  were  calculated  using  the  Ewald  sum  technique 
[17],  which  is  a  well-known  method  for  evaluating  the  electrostatic  interactions  in 
ionic  crystals. 
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The  Ewald  transformation  is  expressed  by  two  convergent  sums.  One  in  real 
space  of  a  short-ranged  potential.  Un, 


r  ,  ^  Z-ZT  erfc<ar"-) 

and  the  other  in  reciprocal  lattice  space  of  periodic  Fourier  domains,  Uh, 


(2a) 


Uh  =■ 


N  N 

ESS 


zizj  q~ 


2-l  .r^rihVo 


exp(-^V)/h2  cos(^Lh.r  ),  (2b) 

q*L'  L 


where  L  is  the  box  length,  a  is  an  arbitrary  parameter,  typically  set  to  5/L,  and 
erfc  is  the  complementary  error  function.  Note  that  h  is  a  reciprocal  wave  vec¬ 
tor  in  units  such  that  its  components  are  integers.  In  the  latter  case  of  reciprocal 
space,  the  order  of  computations  can  be  reduced  by  the  method  suggested  by 
Singer  [18],  in  which  the  double  sum  over  particle  i  and  j  is  simplified  into  two 
single  sums  over  particle  i.  Total  of  618  wave  vectors  were  computed  using  a 
recurrence  relationship  of  complex  arithmetic  to  avoid  repeated  calculations. 


The  system  of  a  hard-core  repulsion  with  continuous  attractive  interactions 
gives  rise  to  problems  for  MD  computer  simulations.  Computational  techniques 
for  trajectory  calculations  are  totally  different  between  the  discontinuous  and  the 
continuous  MD  method.  In  the  hard-core  MD  program  introduced  in  the 
pioneering  work  of  Alder  and  Wainwright  [19],  the  collision  equations  between  all 
possible  colliding  pairs  must  be  solved  before  advancing  their  trajectories.  In 
contrast,  in  the  MD  for  continuous  potentials,  the  time  evolution  of  phase  space 
is  followed  by  solving  the  equations  of  motion  suitably  discretized  in  time. 


Two  distinct  algorithms  can  be  combined  into  the  same  MD  program  by 


returning  to  the  hybrid  "step-by-step"  approach  described  elsewhere  [8,20,21]. 
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The  time  step  interval  is  of  order  of  femto  seconds  which  is  very  small  compared 
to  the  average  time  between  hard-sphere  collisions.  The  first  step  is  identical  to 
the  system  of  conventional  continuous  potentials.  The  system  trajectories  are 
advanced  from  the  current  positions  to  the  next  positions  only  under  the 
influence  of  continuous  forces  without  imposing  the  hard-core  constraints.  In  this 
step,  any  possible  hard-core  collision  occurring  during  a  fixed  time  interval  is 
ignored,  and  this  may  result  in  unrealistic  hard-core  penetration.  The  next  step 
is  then  to  check  whether  or  not  the  pair  distances  are  closer  than  a  hard-sphere 
contact  diameter,  and  the  resulting  configuration  Is  resolved  for  the  overlapping 
pairs.  Because  the  time  interval  is  so  small,  the  majority  of  particles  proceed 
without  colliding  with  each  other.  When  an  overlap  is  found  in  a  time  step,  the 
step  is  repeated,  but  this  time  the  particle  velocities  are  assumed  to  be  constant. 
Under  this  assumption,  the  particles  are  advanced  until  the  overlapping  pair  col¬ 
lides  and  the  algebraic  equations  of  colliding  hard-spheres  are  used  to  find  the 
postcollision  positions  and  velocities.  The  algorithm  then  returns  to  the  continu¬ 
ous  potential  NID  sequence.  This  approach  has  been  applied  for  the  various  sys¬ 
tems  of  model  fluids  and  did  not  show  any  measurable  inaccuracy  for  the  ther¬ 
modynamic  and  transport  properties  of  such  systems  [8,20.21].  A  logic  flow  chart 
for  this  MD  method  is  illustrated  in  Figure  1. 

For  the  computational  efficiency,  it  is  appropriate  to  eliminate  any  redundant 
calculations  and  this  can  be  done  by  the  construction  of  the  collider  table  to 
speed  up  the  searching  routine.  Before  entering  the  force  loop  in  the  first  step, 
the  maximum  velocity  was  scanned  to  decide  the  maximum  cut-off  distance  for 
the  collider  table.  The  table  was  updated  at  each  step  to  ensure  possible  collid¬ 
ing  pairs.  Only  the  pairs  within  this  maximum  distance  were  considered  to 
examine  the  colliding  pairs  in  the  next  step,  rather  than  all  possible  pairs. 
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The  equations  of  motion  were  intergrated  using  the  leapfrog  version  of  Verlet 
algorithm  [2*2],  and  the  velocities  were  scaled  at  each  time  step  to  maintain  con¬ 
stant  temperature  in  the  manner  described  by  Berendsen  et  al.  [23].  In  addition, 
the  initial  configurations  were  generated  by  randomly  inserting  particles  to  assist 
in  the  equilibration  of  the  system.  Configurations  were  aged,  or  equilibrated,  for 
6.000  steps  before  accumulating  data  and  the  resulting  ensemble  averages  were 
obtained  during  the  final  -10.000  steps.  The  MD  algorithm  implemented  here  has 
been  tested  in  a  number  of  ways.  When  the  ionic  charges  were  assigned  to  a 
value  of  zero,  the  results  faithfully  reproduced  the  thermodynamic  and  transport 
properties  of  the  pure  hard-sphere  system  available  in  the  literature  [2-4].  The 
results  obtained  from  the  MD  simulations  for  a  few  selected  runs  of  1-1  electro¬ 
lytes  were  also  compared  to  the  MC  [9]  and  MD  [8]  calculations.  The  good  agree¬ 
ment  with  the  previous  data  again  confirmed  the  quality  of  our  MD  method.  All 
production  calculations  were  performed  on  a  Cray-2  supercomputer  at  Minnesota 
Supercomputer  Center.  Extensive  use  was  made  of  vectorization  and  optimiza¬ 
tion. 


3.  RESULTS  AND  DISCUSSION 


The  thermodynamic  and  transport  properties  of  1-3  electrolyte  solutions 
obtained  from  the  MD  simulations  are  presented  in  Table  1.  The  range  of  con¬ 
centration  investigated  here  was  from  0.02  \ft  to  6  Mt.  where  Mt  is  the  total  ion 
concentration  in  units  of  mole  per  liter,  but  not  the  stoichiometric  concentration, 
Ms.  (Note  that  Mt  =  4  Ms  =  21.6318  nd3,  where  nd3  is  the  total  number  den¬ 
sity.)  The  self-diffusion  coefficients  were  calculated  from  the  integration  of  the 
velocity  autocorrelation  function  using  the  Green-Kubo  relations  and  from  the 
slope  of  the  mean-square  displacement  versus  time  using  the  Einstein  equations. 
The  two  methods  gave  results  in  good  agreement,  typically  less  than  3 Cc 
difference,  and  the  self-diffusion  coefficients  in  Table  1  are  averages  of  the  values 
found  by  the  two  methods.  In  this  table,  we  also  report  the  ion/ion  collision  fre¬ 
quencies  in  columns  7  to  9. 

In  Figures  2  and  3  we  illustrate  the  concentration  dependence  of  the  reduced 
configurational  energy,  U/NkT,  and  the  reduced  osmotic  pressure,  PV/NkT. 
respectively.  Also  shown  in  these  figures  are  the  results  obtained  with  the 
exponential  form  of  Debye-Huckel  theory  (DHX)  [9,15],  the  mean  spherical 
approximation  (MSA)  [3],  and  those  reported  in  Ref.  [4]  for  the  hypernetted 
chain  theory  (HNC).  The  thermodynamic  properties  can  be  expressed  in  terms 
of  pair  distribution  functions,  g(J (r),  between  ionic  species  i  and  j. 

The  configurational  energy  of  the  fluid  is 

=  ■^LEE/U.)(r)g1J(r)r2dr,  (3) 

and,  the  virial  expression  for  the  osmotic  pressure  is 
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PV  =  U  2?rnd 

XkT  +  3NkT  3 

where  X,  is  the~mole  fraction  of  species  i  and 
pair  function  for  species  i  and  j. 


-££X.xlg|J(d).  (4) 

i  j 

gIJ(d)  is  the  contact  value  of  the 


We  first  focus  our  attention  on  the  results  obtained  from  the  DHX  theory. 
This  simple  approximation  is  inadequate  over  most  of  the  range  of  the  concentra¬ 
tions  investigated  in  this  work.  The  discrepancy  is  gradually  amplified  with 
increasing  concentration.  Although  the  poor  results  are  observed  here,  we  point 
out  that  for  the  1-1  electrolyte  system  and  Mt  <  2.0  the  DHX  energy  predictions 
are  in  reasonable  qualitative  agreement  with  previous  simulation  results  [8].  The 
main  failure  of  the  DHX  approximation  is  that  it  misses  the  structure  of  the  pair 
distribution  functions.  The  oscillatory  behavior  occurring  at  higher  concentra¬ 
tions  due  to  the  hard-sphere  exclusion  volume  is  totally  ignored  in  this  theory. 
The  contact  values  of  the  pair  distribution  functions  of  the  DHX  theory  are  also 
inaccurate. 


This  last  point  is  illustrated  in  Figures  4a  and  4b.  Pair  distribution  functions 
were  computed  by  sorting  the  relative  pair  distances  into  the  equal  radial  incre¬ 
ments  of  width  O.Old.  Even  at  the  low  concentration  of  Mt  =  0.2  (Figure  4a), 
slightly  more  H —  and  +4-  ion  pairings  at  the  short  range  of  distances  were  found 
in  the  MD  result  compared  with  the  DHX  pair  functions.  The  sharp  peak  near 
the  contact  pjuvi  in  g+_  indicates  the  strong  tendency  to  the  formation  of 
unfike-pairs.  As  seen  in  Figure  4b  (Mt  =  4.0),  however,  the  DHX  model 

significantly  overestimates  --  pairings  in  g while  the  underestimation  of  the 

Coulombic  attraction  between  unlike-pairs  results  in  much  lower  value  in  g+_. 
In  the  MC  calculations  for  the  2-2  electrolytes  [9,11],  the  presence  of  linear  ion 
triplets  was  manifested  by  a  local  maximum  in  like-pair  functions  near  r  =  2d. 
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Similar  results  were  observed,  but.  in  the  ease  of  our  1-3  system,  the  maximum 
peak  in  g__  was  shifted  to  approximately  r  =  2.5d.  The  salient  feature 

displayed  here  is’ a  noticeable  trend  to  the  formation  of  larger  ionic  clusterings 
for  higher  charged  systems. 

The  MSA  and  HNC  theories  are  both  based  on  the  Ornstein-Zernike  integral 
equation  along  with  specific  closure  relationships.  The  MSA  model  has  a  great 
advantage;  it  is  the  only  theory  that  can  be  solved  analytically  for  the  system  of 
charged  hard-spheres  [3].  .\s  an  extension  to  the  Debve-Huckel  theory,  in  the 
MSA  theory,  the  hard-sphere  oscillations  due  to  the  ion-cores  have  been  taken 
into  account.  In  the  limit  of  point  charges,  the  MSA  becomes  identical  to  the 
Debye-Huckel  theory.  .As  illustrated  in  Figure  2  and  3,  the  MSA  results  for  the 
energy  and  for  the  osmotic  pressure  are  in  relatively  good  agreement  with  the 
MD  data.  However,  the  MSA  pair  distribution  functions  are  known  to  be  poor. 
Particularly  for  a  range  of  low  concentration,  the  linearized  nature  of  the  MSA 
model  leads  to  the  unphysical  negative  values  among  the  like-pair  functions  as  in 
the  Debye-Huckel  theory.  The  MSA  pair  functions  are  only  symmetrically 
disposed  around  those  for  the  hard-sphere  pair  function.  The  better  agreement 
for  pressure  than  for  energy  in  the  DHX  model  seems  to  be  coincidental  and 
caused  by  a  cancellation  of  errors  between  energy  and  collision  contributions  (the 
second  and  the  third  term  in  the  virial  pressure  in  eq.  (4),  respectively).  For  Mt 
=  2.0,  for  instance,  an  underestimation  in  g+_(d)  by  about  30%  is  almost 
equally  balanced  by  an  overestimate  in  energy.  The  HNC  results  for  the  system 
of  1-3  electrolytes  [4]  predict  the  configurational  energy  and  the  osmotic  pressure 
in  very  good  agreement  to  the  MD  results  within  statistical  errors.  Although  spe¬ 
cial  iterative  techniques  are  necessary  in  order  to  reach  convergence  [4],  HNC 
theory  has  also  shown  to  be  successful  for  the  1:1  and  2:2  electrolytes  [9,11]. 
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The  velocity  autocorrelation  functions  (VACF)  can  provide  useful  insights 
into  ion  dynamics  and  transport.  The  normalized  VACF  is  defined  as  a  function 
of  time.  t.  by 


VACF  =  ^  £0,(0)  .  v,(t)>  /  <v,-(0)>.  (5) 

i-  l 

where  the  symbol  <  >  denotes  the  average  over  an  equilibrium  ensemble. 
Those  time  correlation  functions  are.  in  principle,  measurable  quantities  from 
inelastic  scattering  experiments. 

In  Figures  5a  to  5c,  we  plot  the  normalized  VACFs  for  a  few  selected  runs  to 
illustrate  the  manner  in  which  those  functions  change  with  increasing  concentra¬ 
tion.  The  solid  curves  correspond  to  the  VACFs  for  positively  charged  ions, 
VACF(+),  and,  respectively,  the  dotted  curves,  for  negatively  charged  ions, 
VACF(— ).  The  VACFs  for  lower  concentration  exhibit  a  stronger  velocity  auto¬ 
correlation  than  those  for  higher  concentration.  The  primary  mechanism  of 
decay  of  the  time  correlation  functions  is  the  hard-sphere  collision,  in  which  col¬ 
liding  particles  rapidly  forget  their  initial  velocities  through  successive  collision. 
In  the  high  concentration  regime,  the  hard-core  repulsive  collisions  are  expected 
to  play  a  dominant  role  in  determining  the  dynamical  properties  of  the  system. 
In  agreement  with  this,  the  resulting  VACF(+)s  exhibit  the  exponential  behavior 
which  characterizes  the  dynamics  of  hard-sphere  fluids  [24]. 

However,  as  shown  by  the  VACF(— )s  in  Figure  5,  there  are  substantial 
differences  between  the  motion  of  the  positive  charges  and  the  highly  charged 
negative  ions.  The  VACF(— )s  decay  more  rapidly  than  the  VACF(+)s.  The 
VACF(— )s  exhibit  non-exponential  behavior.  At  low  or  intermediate  concentra¬ 
tion,  the  long-ranged  Coulombic  potential  increasingly  influences  the  collective 
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motion  of  ionic  fluids  by  the  acceleration  or  retardation  of  velocities  between  col¬ 
liding  pairs.  The  electrostatic  force  enhances  correlations  between  unlike-pairs. 
and  decorrelation  between  like-pairs. 

The  most  striking  feature  of  the  VACF(— )s  is  the  peculiar  peak  displayed  in 
the  range  of  intermediate  concentration.  The  peak  occurs  at  relatively  the  short 
time  delay  between  5  ps  and  10  ps.  At  the  very  low  concentration  regime,  the 
YACF(— )  decays  monotonically.  However,  a  weakly  oscillatory  behavior  starts 
to  emerge  in  the  VACF(— )  for  Mt=0.5  (Figure  5a),  and  is  clearly  apparent  for 
Mt=  1.0  (Figure  5b).  This  can  be  explained  in  terms  of  the  coherent  motion  of  a 
highly  charged  negative  ion,  and  the  positive  ions  surrounding  it.  The  coherent 
motion  in  this  range  of  concentration  results  in  an  increased  persistence  of  velo¬ 
city  of  the  negative  ions.  Because  the  cluster  tends  to  move  coherently  over  a 
period  of  time  longer  than  the  mean  collision  time  between  ions  in  the  cluster, 
this  collective  motion  conducts  the  negative  ion  in  the  direction  of  the  net  drift 
of  the  entire  cluster.  In  the  cases  presented  in  Figures  5a  and  5b,  the  bump  in 
the  VACF(— )  is  located  at  nearly  twice  the  mean  collision  time  between  negative 
and  positive  particles.  This  fact  suggests  that  the  mechanism  responsible  for  this 
velocity  persistence  is  a  double  collision.  First,  a  collision  occurs  between  the 
negative  central  ion  and  the  positive  ions  in  front  of  it.  Then,  a  second  collision 
occurs  when  the  same  negative  ion  is  reached  by  the  positive  ions  which  the  elec¬ 
trostatic  attraction  drag  behind  it.  In  contrast,  as  seen  in  Figure  5c,  these  bumps 
are  reduced  at  the  higher  concentration  of  Mt=  6.0.  In  this  case,  there  is  perhaps 
not  enough  free  volume  for  the  ionic  clusters  to  move  coherently  over  appreciable 
distances. 

In  the  previous  MD  simulations  (8)  for  the  1-1  electrolytes,  the  modified 
Enskog  theory  of  hard-sphere  was  shown  to  predict  the  self-diffusion  coefficients 


reasonably  accurately.  It  is  also  interesting  to  note  that  the  diffusion  constants  of 
positive  ions  for  the  1-3  electrolytes  here  (Table  1)  are  close  to  those  for  the  1-1 
electrolytes  except  for  the  low  concentration  limit  of  Mt  <  0.2.  However,  the 
self-diffusion  coefficients  of  negative  ions  are  smaller  by  a  factor  of  two  or  three 
than  those  of  positive  ions.  An  interpretation  of  this  observation  is  that  the  free 
motion  of  negative  ions  is  likely  to  be  restricted  by  the  formation  of  ionic  aggre¬ 
gates. 

For  a  hard-sphere  system  [26],  the  collision  frequencies,  u,\. ,  are  expressed  in 
terms  of  the  contact  values  of  pair  distribution  functions, 

=  4ihd2  (7rkT/m)1/2gjj(d),  (6) 

where  is  the  number  of  collisions  per  particle  of  component  j,  per  unit  time 
between  particles  of  component  i  and  j.  The  collision  frequencies  determined 
from  the  MD  simulation  are  plotted  in  Figure  6.  For  the  purpose  of  comparison 
with  the  corresponding  hard-sphere  system,  the  predictions  using  the  MD  contact 
values  in  eq.  (6)  are  also  shown  in  this  figure.  The  MD  results  are  seen  to  be  in 
excellent  agreement  with  theory.  This  suggests  that  the  microscopical  behavior 
of  the  primitive  electrolytes  is  much  like  to  the  hard-sphere  dynamics.  In  this 
sense,  the  transport  properties  of  the  primitive  model  electrolytes  are,  at  least 
qualitatively,  related  to  those  for  the  system  of  hard-sphere  fluids. 
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4.  Conclusion 

In  the  present  paper  we  have  reported  the  \H)  simulations  for  the  thermo¬ 
dynamic  and  transport  properties  of  1-3  electrolyte  solutions.  The  results 
obtained  from  computer  simulations  have  been  used  to  assess  the  applicabilities 
of  theoretical  equations,  namely,  the  exponential  form  of  Debye-Huckel  theory, 
the  mean  spherical  approximation,  and  the  hypernetted  chain  theory.  The  HXC 
theory  is  shown  to  be  the  best  approximation  compared  with  the  simulation 
results  both  for  energy  and  for  osmotic  pressure.  For  asymmetrically  charged  1-3 
bulk  electrolytes,  the  presence  of  ionic  clustering  was  observed  in  pair  distribu¬ 
tion  functions.  The  local  maximum  in  g__  was  shifted  to  r  =  2.5d. 

The  MD  results  indicate  that  the  diffusion  processes  are  strongly  dependent 
on  concentration  and  the  solvent  have  a  negligible  effect  on  the  dynamics  of 
solute  ions.  However,  experimental  diffusion  coefficients  in  aqueous  solutions  [27] 
are  nearly  concentration  independent.  Under  certain  conditions  the  primitive 
model  is  known  to  be  a  realistic  description  for  molten  salts.  One  important 
aspect  neglected  in  the  primitive  model  electrolytes  is  the  structural  changes  of 
the  solvent  according  to  the  local  concentration  of  solute  ions.  This  simple  model 
is  an  inadequate  to  describe  the  ionic  structures  and  dynamics  in  aqueous  electro¬ 
lyte  solutions.  With  the  complicated  vibrational  and  orientational  motions  of 
water  molecules,  numerous  factors  are  involved  to  ascertain  the  principal  features 
governing  the  transport  properties  of  aqueous  solutions.  Consequently,  the  satis¬ 
factory  model  potentials  describing  the  detailed  ionic  motions  in  aqueous  solu¬ 
tions  are  required  for  further  work. 
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Table  1.  System  characteristics  and  MD  results 
for  1-3  Primitive  Model  Electrolytes 


Mt 

-U/NkT 

PV/N'kT 

D+ 

D_ 

u'  +  + 

W _ 

W-  + 

(moles//) 

(10  4cm2/s) 

(10-4cm2/s) 

( 10 1  °s—  1 ) 

(1010s-‘> 

(10I0s-1) 

0.02 

0.5169  (0.0437)' 

0.8558 

315.21 

113.84 

0.068 

0.000 

1.929 

0.2 

1.1983  (0.0543) 

0.7433 

57  42 

22  26 

1.042 

0  000 

11.119 

0.5 

1.5010  (0.0495) 

0.7239 

28.54 

12.94 

2.801 

0.000 

17.185 

1. 

1.7423  (0.0495) 

0.7370 

17.24 

8.817 

5.977 

0.007 

23.603 

o 

2.0022  (0.0469) 

0.8110 

9.574 

5.553 

13.142 

0.014 

33.174 

4. 

2.2978  (0.0438) 

1.0293 

5.151 

3.288 

31.182 

0.019 

50.895 

6. 

2.4957  (0.0405) 

1.3495 

3.490 

1.699 

55.679 

0.167 

70.319 

*  The  values  in  parentheses  indicate  uncertainties  in  the  MD  simulations. 
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Figure  Captions 

1.  Schematic  of  the  niolecular  dynamics  procedure. 

2.  Reduced  configurational  energy  of  1-3  primitive  model  electrolytes  as  a  function  of  concentra¬ 
tion.  □  NED  results;  - DHX: - MSA;  -  •  -  HNC[4]. 

3.  Reduced  osmotic  pressure  of  1-3  primitive  model  electrolytes  as  a  function  of  concentration. 
The  symbols  have  the  same  meaning  as  in  Figure  2. 

4.  Pair  distribution  functions  for  like-pairs  and  unlike-pairs.  (a)  MT=  0.2;  (b)  Mt=  4.0.  Also 
shown  as  dotted  curves  are  DHX  predictions. 

5.  Normalized  velocity  autocorrelation  function  vs.  time  at  three  different  concentrations:  a) 
Mt=  0.5,  b)  1.0,  and  c)  6.0.  The  solid  curves  correspond  to  the  velocity  autocorrelation  function 
of  positive  ions  (VACF(+ ))  and  dotted  curves  to  that  of  negative  ions  (VACF(— )). 

6.  Semilogarithmic  plot  of  collision  frequency  as  a  function  of  concentration. 

O  MD  results  for  k/+  +  ;  ■  MD  results  for  iU _ The  dotted  curves  are  predictions  by  Eq.  (6) 

using  the  MD  contact  values  of  pair  distribution  functions. 
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